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Abstract 



of rather general (3+l)-dimensional data sets can be computed in order 
to solve the Einstein equations. After a general introduction, three topics 
of current interest are briefly reviewed: binary black hole mergers, the 
evolution of strong gravitational waves, and shift conditions for neutron 
star binaries. 
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$h ! 1 Introduction 

W). 

In this short review I describe work concerned with one of the central issues of 
numerical relativity, the solution of the two body evolution problem of general 
relativity. After a short introduction to (3+l)-dimensional numerical relativity, 
I briefly discuss recent progress on binary black hole mergers, the evolution of 
strong gravitational waves, and shift conditions for neutron star binaries. 

As opposed to Newtonian theory, where the Kepler ellipses provide an astro- 
physically relevant example for the analytic solution of the two body problem, 
in Einsteinian gravity there are no corresponding exact solutions. The failure 
of Einstein's theory to lead to stable orbits is due to the fact that in general 
two orbiting bodies will emit gravitational waves that carry away energy and 
momentum from the system, leading to an inspiral. But, of course, this "leak" 
is not considered to be detrimental. Gravitational waves are one of the most 
interesting new phenomena introduced by general relativity that will open a 
new window into the universe through gravitational wave astronomy, e.g. jjj]. 

The evolution of a two body gravitational system, for example a binary 
black hole system (which can be constructed as a vacuum system and avoids 
additional complication due to matter sources), can be divided into at least 
three phases. For sufficiently large separation of the two black holes there is a 
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slow inspiral phase with many orbits, followed by a very brief violent merger 
phase that leads to a single, distorted black hole that after a short ring-down 
phase settles down to a final stationary black hole. For the initial and final phase 
rather well understood approximation schemes are available, i.e. post-newtonian 
calculations for the slow inspiral of two point masses (e.g. [|| [|) and the close 
limit approximation for the ring-down of a single distorted black hole(e.g. 0)). 
For a full treatment of the strongly non-linear, fully general rclativistic phase one 
has to turn to computer simulations to obtain (again approximate) numerical 
answers. 

Each phase leads to a characteristic gravitational wave signal. At this time 
several gravitational wave detectors are being built world-wide that should for 
the first time make the direct measurement of gravitational waves possible. The 
prediction and analysis of future signals is the main motivation for studies of 
binary systems in numerical relativity. While certainly the primary motivation, 
let me add that even if it had no direct, in the near future measurable observa- 
tional consequences, we should solve a basic problem like the two body problem 
of general relativity. 

The article is organized as follows. In Sec. |[ a brief history of black hole 
evolutions in numerical relativity in 2+1 (axisymmetry) and 3+1 dimensions 
is given. In Sec. |^, the evolution problem of numerical relativity is introduced 
in its 3+1 form, leading to three main issues: initial data, evolution, analysis. 
Initial data is computed on a three-dimensional hypersurface, which is evolved 
in time, and at various times analysis like identifying the black hole horizons 
and gravitational wave extraction is carried out. The coordinate problem of 
numerical relativity is emphasized with the choice of slicing function, the lapse, 
as an example. The skeleton of a typical black hole evolution is discussed. In Sec. 
|], the "Cactus" code, a computational infrastructure for numerical relativity 
and relativistic astrophysics, is described. 

After this general introduction, we discuss several examples for recent pro- 
gress in (3+l)-dimensional numerical relativity. In Sec. ^, current (summer of 
1999) binary black hole simulations are presented. The holes start out close to 
each other and evolve through a plunge rather than an orbit (a grazing collision). 
Achievable evolution time is now about 30M (M the mass of the final black 
hole), which for the first time allows the extraction of wave forms. In Sec. [| 
a discussion of strong wave evolutions is included because strong waves play a 
role in black hole mergers and these studies provided the proving ground for a 
new evolution scheme discussed in Sec. 3.1. In Sec. [t], the minimal distortion 
shift condition is described. Lapse and shift specify the coordinate gauge, and 
in all simulations mentioned so far the shift has been zero, but for systems 
with rotation a shift condition will be essential. The example presented is 
minimal distortion shift for a binary neutron star system, which for this purpose 
is simpler than black holes because there are no special inner boundaries. 

Sec. H concludes this brief review, pointing out again those issues and tech- 
niques that will be important for the numerical simulation of binary black holes 
for several orbits lasting for 100-1000M with results that are relevant for grav- 
itational wave astronomy. 
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2 History of black hole simulations in numerical 
relativity 

In this section I endeavor to give a necessarily very short but in its highlights 
complete exposition of the literature on numerical black hole evolutions, con- 
centrating mostly on work that implements the complete black hole evolution 
problem (data, evolution, analysis) for the full Einstein equations in vacuum. 
Matter occurs in a few places but only as a means to form black holes. Clearly, 
there is a large and important body of work concerned with all the different, 
separate aspects of and methods for black hole evolutions as outlined in Sec. |[ 
Still, this allows us to sketch the history of the field. 



2.1 2+1 dimensions 

After some early attempts (5) , it was the work by Smarr 0, fjj and Eppley || on 
the head-on collision of two equal mass Misner black holes |i| which basically 
founded the field of numerical relativity as a subject of computational physics. 
Axisymmetric head-on collisions allow significant savings in computational cost 
when formulated with two spatial and one time coordinate (2+1 dimensions), 
although this excludes the possibility of orbiting black holes and radiation of 
angular momentum. Many of the key techniques that are still in use today stem 
from that period of the sixties and seventies (see (To) for the definitive review). 

The beginning of the nineties saw a surge of activity when more power- 
ful computers, improved codes and methods allowed significant advances. The 
axisymmetric collision of black holes, either formed by particles Jll| or imple- 
mented as Misner data [l2| [l3| , was repeated complete with horizon finding and 
wave extraction. It is remarkable how the crude results for the wave emission 
of |(| [?| where confirmed in |Q. Another highlight is certainly the numerical 
computation of the "pair of pants" picture for a black hole merger |l4j , which 
was a result of the US Binary Black Hole Grand Challenge Alliance. Head-on 
collisions in axisymmetry continue to improve, see the recent work on unequal 
mass configurations . 

Another interesting system in axisymmetry are rotating black holes. In ]l7t , 
particles collapse to form a black hole with rotation and toroidal event horizon. 
In |l^| , a Kerr black hole distorted by a gravitational wave is evolved. Matter 
plus rotating black hole systems are also studied in |20| . 

A traditional topic in numerical black hole studies is that of black hole for- 
mation, of which I want to mention only the following recent references that are 
of relevance to this article. The formation of naked singularities was examined 
in Furthermore, in [^2| the collapse of gravitational waves to a black hole 
is demonstrated. A surprise was that even (l+l)-dimensional, spherically sym- 
metric black hole systems are far from being trivial, as the rich set of critical 
phenomena discovered by Choptuik J|3| showed (see e.g. Q for a review). In 
2+1 dimensions, the only critical collapse studies so far are those of [ p5| . 



3 



2.2 3+1 dimensions 



Numerical relativity of black holes in 3+1 dimensions was initiated in 1995 with 
evolutions of a Schwarzschild black hole with singularity avoiding slicing on a 
Cartesian grid pq] . Achieved run time is about 30 M. At the same time the 
first (3+l)-dimensional wave simulations were carried out p8| . In Sec. |[ I 
comment on the collapse of non-axisymmetric waves to a black hole [g9| . 

Returning to our main topic, the evolution of a Schwarzschild black hole 
with a non- vanishing shift vector was studied in ]3C|], compare Sec. @. In @, 
adaptive mesh refinement techniques, made famous in numerical relativity by 
[ p3| , were applied for the first time to 3+1 relativity, also for a Schwarzschild 
black hole. By now, evolutions for the Schwarzschild spacetime are standard 
code test, e.g. [|32[ . Further studies of single black holes include the distorted 
black holes in J33[ |34], |35[ , which provided the first detailed tests of wave extrac- 
tion in 3+1 dimensions. The Black Hole Grand Challenge Alliance performed 
the longest, stable evolution of a single black hole so far, reaching about 60M 
for a standard Cauchy evolution with black hole excision and a boosted black 
hole |56|, and essentially achieved complete stability (> 60000M) with a char- 
acteristic evolution code, which is tailored to the one black hole problem but 
can also treat small distortions, and for the first time a black hole that moves 
across the grid © H !§• 

Binary black hole evolutions are pushing the limits of what is currently 
possible. Some results for the evolution of the axisymmetric Misner data set with 
the 3+1 code of [^6| with singularity avoiding slicing are reported in The 
first true (3+l)-dimensional binary black hole evolutions, the grazing collision 



of nearby spinning and moving black holes, was performed in |41 . This sets the 
stage for the recent binary black hole simulations of Sec. ||, but first we want to 
discuss some of the basic issues in numerical relativity. 



3 Anatomy of a numerical relativity simulation 
3.1 3+1 formulation 

The Arnowitt-Deser- Misner (ADM) equations |^2|, |l(J are one of the possibilities 
to rewrite the Einstein equations as an initial value problem for spatial hyper- 
surfaces. The dynamical fields of the ADM formulation are a 3-metric g ab and 
its extrinsic curvature K ab on a 3-manifold S, both depending on space (points 
in S) and a time parameter, t. The foliation of the 4-dimensional spacetime 
into hypersurfaces £ is characterized in the usual way by a lapse function a and 
a shift vector (3 a . The Einstein equations for vacuum become 

{d t -Cp)g a b = -2aK ab (1) 

{dt-£ )K ab = -D a D b a + a(R ab -2K ac K c b + K ab K) (2) 

= D b {K ab - g ab K)=V ai (3) 

= R - K ab K ab + K 2 = H, (4) 
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where R a b is the 3-Ricci tensor, R the Ricci scalar, K the trace of the extrinsic 
curvature, Cp the Lie derivative for (3, and D a the covariant derivative compati- 
ble with the 3-metric. One obtains evolution equations for the metric variables, 
(Q) and (^) , and constraint equations that do not contain time derivatives of g a b 
or K a b, the momentum constraint (|3|) and the Hamiltonian constraint (Q). 

These equations are well known, but displaying them explicitly allows me 
to make a number of basic observations. First of all, these are comparatively 
simple equations. Even though writing out all the terms in the index contrac- 
tions, in the definition of R a b and the covariant derivative leads to on the order 
of 1000 floating point operations per point for a typical finite difference repre- 
sentation of (jl|) and (||), this can be easily dealt with computationally and is 
not one of the fundamental problems of black hole evolutions. Still, a numerical 
implementation requires some thought and hard work, see Sec. ^. 

Notice that lapse and shift appear in the evolution equations (as of course 
they have to) and have to be specified as part of the evolution problem. Choosing 
lapse and shift fixes the coordinate gauge for the evolutions, and is one of the 



key problems for numerical evolutions, see Sec. 3.2 



The constraint equations imply that specifying initial data g a b and K a b on 
a hypersurface £ involves in general solving the constraints numerically. If 
the constraints are satisfied initially, they will remain satisfied for a well-posed 
evolution system, but this is an analytic statement that is only approximately 
true numerically 

Finally, the ADM equations do not define a hyperbolic evolution system 
(e.g. Q]), and it is not clear to what extent the original ADM equations can 
lead to a numerically stable evolution system. The issue of stability has to be 
addressed on two levels. First, a well-posed evolution system is one for which 
existence, uniqueness and stability of a solution for at least finite time intervals 
can be shown, which is true for example for hyperbolic systems. However, 
in general stability does not rule out exponentially growing modes (this may 
be the solution one is looking for). Second, the numerical implentation of an 
analytically stable system does not trivially lead to stable numerical evolutions 
(e.g. the finite-differenced equations may have exponentially growing solutions 
which are not present analytically). Also note that important stability issues 
arise at the boundaries of the computational domain. 

Finding stable evolution systems is perhaps the other key issue in numerical 
relativity besides the choice of coordinate problem. For an excellent review 
of first order hyperbolic systems for relativity see j44|, but other systems are 
of interest, too. One of the important developments of the last year was the 
demonstration by Baumgarte and Shapiro that a conformal, trace-split 
version of the ADM system very much like the system used by Shibata and 
Nakumara in p7| ], is significantly more stable (numerically) than the ADM 
equations for weak fields and some algebraic slicings. This BSSN system can 
also be understood as a second order, conformal version of the Bona-Masso 
system BO] . First order hyperbolic versions were given in p7L H , although the 
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Figure 1: Schwarzschild black hole in Novikov coordinates, M = 1. 



BSSN system as it stands is not hyperbolic. The BSSN variables are 



= ln(det 5 )/12, (5) 

K = g ab K ab , (6) 

~g ab = e-^gab, (7) 

A ab = e~^{K ab -g a bK/i), (8) 

f c = f c ab g ab , (9) 



so that Aetg = 1 and trA a b = 0. Furthermore, introducing T c leads on the 
right-hand-side of equation (^) to an elliptic expression in derivatives of g a b, 
i.e. the corresponding BSSN equation has the character of a wave equation. 
However, including the evolution equation for r c appears to spoil hyperbolicity. 
Nevertheless, the BSSN system has very nice stability properties, and some 
suggestions about why this may be the case are made in (49|. Several BSSN 
evolutions have now been reported, for strong waves and maximal slicing |p9~f 
(Sec. ^) and also for matter evolutions 

3.2 Schwarzschild as an example for a typical black hole 
evolution problem 

Moving on to the prototypical example for a black hole evolution, consider Fig. 
[j], which shows the Schwarzschild spacetime for a static, spherically symmetric 
black hole in Novikov coordinates |33|, The coordinates are chosen such 
that freely falling observers that start at rest at time r = follow constant 
R* lines. The Schwarzschild radius r is related to R* at r = through R* = 
(r/(2M) — l) 1 / 2 . Several constant r lines are shown: the physical singularity 
at r = 0, the event horizon at r — 2M, and note how the lines r — AM and 
r = 6M curve outwards which corresponds to the radial infall of the observers 
with constant R* . 
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To set up an evolution problem we can choose the slice r = as initial 
hypersurface with g a b and K a b derived from the Schwarzschild four-metric {g a b 
and K a b therefore solve the constraints). Note that the physical singularity 
is to the future of this slice and does not show in g a b and K a b. To perform 
an evolution, we have to specify lapse and shift. The Novikov coordinates 
correspond to geodesic slicing, a = 1, and vanishing shift, (3 a = 0. Concretely, 
consider a numerical grid at r = extending from R* = to R* = 2 1 / 2 . In the 
first quadrant of the figure we have shown how this initial slice moves through 
the spacetime for geodesic slicing with vanishing shift. Without precaution a 
numerical code will crash at r = irM when the point at R* = reaches the 
singularity (this "crash test" has in fact been used as a first crude code test 
©0). 

As shown in the figure, one can imagine evolving beyond r = nM by cutting 
out from the slice what is inside the event horizon of the black hole, which will 
not affect the outside of the black hole anyway. "Black hole excision" techniques 
[ |55[ |5(| are a very promising approach to black hole evolutions, although in 3+1 
dimensions there remain certain stability problems to be resolved for binary 
black holes. Black hole excision usually involves a non-trivial choice of lapse 
and shift. 

With black hole excision not quite ready yet, it is so-called singularity avoid- 
ing slicings that have been most widely used. Assume (3 a = 0. Primary examples 
are maximal slicing (K = initially and Aq = aK a bK ab so that dtK = 0), and 
so-called "1+log" slicings (e.g. d t a = —aK/2). In the second quadrant of Fig. 
|l|, we show a typical example (hand-drawn, while the rest of the figure was com- 
puted). At the center, evolution slows down, while for large radii the evolution 
marches on with a — 1 at infinity. Obviously, a numerical problem will occur 
in between, which is referred to as grid-stretching, and which is reflected in 
growing sharp peaks in the radial-radial component of the metric. Singularity 
avoiding slicings have this fatal problem built in, but they allow us to compute 
evolutions up to 30M — 100M, which is barely sufficient for certain black hole 
collision and ring-down wave forms. 

For completeness let me also mention the possibility of using hyperboloidal 
slices, or null-slices, or null-slices matched to the spatial slices, which when 
applicable cover more of the interesting space time in the wave region, e.g. 
[0 HI' H' HI 13' H' IE @- Characteristic matching is also useful near an 
excision boundary. 

The key point to note is that the choice of coordinates in relativity is a 
more fundamental problem than, say, the choice of spherical over Cartesian 
coordinates for computational convenience. There is no simple canonical choice 
like a = 1 and (3 a = that works in reasonably general situations. Even if there 
are no black holes, geodesic slicing fails due to geodesic focusing. It appears to 
be the case that one has to determine lapse and shift dynamically during the 
evolution by some geometric principle as functions of the metric and extrinsic 
curvature. 
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3.3 Anatomy of a black hole simulation 



After discussing 3+1 formulations in general and a specific black hole example, 
let us list the components of a numerical relativity evolution, with the binary 
black hole problem as example. 

3.3.1 Initial Data 

• Choice of hypersurface 

Simplest choice is R 3 for non-black hole data. Black hole data can be, for 
example, of Misner type with an isometry boundary condition at spheres 
representing the throats of the holes, R 3 — sphere s f6q| , or Brill-Lindquist 
type data based on a punctured R 3 , R 3 — points |6(f| . 

• Solution to constraints 

There are four constraint equations that restrict the choice of 12 com- 
ponents in g a b and K a b- The most common approach is the conformal 
method ©|| 0. 

3.3.2 Evolution 

• Variables and evolution system 

There are many different choices that can be roughly divided into ADM 
like systems that are of second order [^7], [l|| , and first order systems that 
are constructed to obtain hyperbolicity, e.g. Ji4| . 

• Choice of coordinates (gauge choice) 

Typically a vanishing shift is used, but see Sec. [?]. For the lapse, as 
explained above, algebraic and elliptic conditions are in use. 

• Physical singularities 

Smooth, regular initial data may develop physical singularities which are 
features of black hole spacetimes. Physical singularities can be avoided by 
choice of slicing, or removed from the grid through black hole excision. 

• Coordinate singularities 

Dynamical determination of lapse and shift may lead to coordinate patho- 
logies, in particular for algebraic slicings (e.g. (6^, [70| . Elliptic conditions 
are sometimes preferable, although they are computationally much more 
expensive. 

• Outer boundary condition 

Asymptotic flatness can be assumed, which implies fall-off conditions for 
the fields. For run times on the order of 100M for a typical singularity 
avoiding Cauchy evolution, a radiative boundary condition is sufficiently 
accurate and stable, e.g. For a more sophisticated scheme see |f7lf . 

Also there are two well developed approaches in which the numerical grid 
does not end at a finite radius but extends to future null infinity, either 
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by matching to a characteristic code at finite radiu s J62| |63|, o r smoothly 
without matching via a conformal transformation |57|, |58|, f59j [72j. (Note 
that in 3+1 dimensions it is no longer straightforward to use a logarithmic 
radial coordinate as is conventionally done in axisymmetry.) 

Inner boundary conditions 

As discussed above, black hole excision leads to a particular inner bound- 
ary. As for initial data construction, the inner boundary for slices in a 
black hole spacetime may be spherical or point-like. For short term evo- 
lutions, the numerical slice can cover the inner asymptotically flat regions 
of the black holes if the resulting coordinate singularities are treated with 
the puncture method for evolution |2q, 41 



3.3.3 Analysis 

• Tensor components 

The raw output of a computer code will be the components of its basic 
variables, e.g. g a b, K a b, a, and (3 a , all other information is computed from 
these. Interesting local quantities include Riemann curvature invariants / 
and J and the Newman-Penrose invariants ?po through i\>^. 

• Black hole horizons 

The event horizon of black holes is a spacetime concept and can be found 
approximately if a sufficiently large spacetime slab has been computed, 
e.g. fbfll . The apparent horizon is a notion intrinsic to the hypersurface 
(and is therefore slicing dependent). It is defined as the union of outer- 
most marginally trapped surfaces, i.e. surfaces for which the expansion 
of outgoing null rays vanishes. Trapped surfaces are linked to the exis- 
tence of black holes through the singularity theorems. See e.g. j73| and 
references therein for numerical issues. 

• Wave extraction 

Wave forms can be computed reliably at finite but large radius using the 
first order gauge invariant approach of |74| , as recently demonstrated for 
3+1 dimensions in J35|. In approaches that make future null infinity part 
of the numerical grid, see above, wave extraction is much more direct. 



4 Implementation of a numerical relativity sim- 
ulation 

As should be evident from the previous section, numerical relativity poses a 
complex scientific problem that translates into a challenging software engineer- 
ing problem. Here I want to discuss "Cactus", a code that is developed and 
used at the Albert-Einstein-Institut (AEI, the Max-Planck-Institut fur Gravi- 
tationsphysik), and several other institutions |7f| f7(], 77 . 
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Referring to Jr5| , the cactus code is a freely available modular portable 
and manageable environment for collaboratively developing high-performance 
multidimensional numerical simulations. Cactus provides a powerful applica- 
tion programming interface based on user modules (thorns) that plug into a 
compact core (flesh). Cactus is composed of modules that are independent of 
relativity, and of modules designed for relativity. The Cactus Computational 
Tool Kit supports a variety of supercomputing architectures and clusters, im- 
plements MPI-bascd parallelism for finite difference grids, several input/output 
layers, elliptic solvers, metacomputing, distributed computing, and visualiza- 
tion tools. Fixed and adaptive mesh refinement is under development. Cactus 
significantly enhances collaborative development by providing code sharing via 
CVS and defining appropriate interfaces for code combination. A large number 
of physics modules or thorns are available for numerical relativity and astro- 
physical applications, e.g. there are thorns for initial data, evolution routines, 
and data analysis. The first version of Cactus was created by J. Masso and P. 
Walker, and has been available for testing since April, 1997 {78j. The Cactus 
Computational Tool Kit (Cactus 4.0) saw its first public release as a community 
code in July, 1999. It is actively supported by a cactus maintenance team, and 
there is good documentation. Cactus is a "third generation" code, going back 
to the "G" |26| and "H" |2^] codes. The key step taken forward is the massive 
investment in the collaborative infrastructure, which is now beginning to pay 
off. For many more details, see |75|, f76[ . 

So what does all of this mean? Suppose you want to run a black hole 
simulation. Cactus is not a high-level science tool where you get an executable 
with graphical user interface to input, say, the black hole masses and off it goes. 
Numerical relativity is still too experimental for that. At its heart, Cactus 
is a large collection of source files together with a sophisticated make system. 
The user decides what sources to include, then compiles the code. Runs are 
controlled by a text file containing parameters, e.g. for the grid size and the 
black hole masses. The sophistication lies in the ease with which code can be 
changed or added by single users without affecting functionality provided by 
others. Suppose a users wants to add a routine that computes the determinant 
of the metric. A new thorn is created with the source code, e.g. 20 lines of C 
or Fortran, and with files that inform Cactus about new parameters, new grid 
functions (say an array of reals with the name "detg"), and tell cactus when to 
call the new routine (in this case whenever analysis is done). This is work to 
be done by the thorn writer, but he or she gets the rest for free: set-up of a 
numerical grid, storage for g a b and its determinant, evolution of g a b according 
to, say, the ADM equations, parallel execution, input, output, adaptive mesh 
refinement, etc. When submitted to the Cactus code repository, any Cactus 
user can now make use of "detg" . 

Taking the viewpoint of a physicist, the Cactus infrastructure takes care of 
many computer tasks that often distract from science. To say that a simula- 
tion was carried out with Cactus can refer to Cactus, the Computational Tool 
Kit, in the same way that credit is given to MPI for parallelism, or Mathc- 
matica or Maple for symbolic computation. Cactus is successful if the science 
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outweighs the infrastructure. In order that Cactus does not remain a faceless 
collection of source code, I would like to give several science examples and also 
to mention at least a few names in connection with science projects. Work 
on hyperbolic methods in numerical relativity |79| was done by J. Masso, P. 
Walker, and others. A project on black hole excision techniques has been im- 
plemented as "Agave" by S. Brandt, M. Huq, P. Laguna, and others at Penn 
State University, which uses Cactus mainly for parallelism. Furthermore, the 
NASA Neutron Star Grand Challenge Project of W.-M. Suen, E. Seidel, and 
others, develops the so-called GR3D code f80fl , which is a version of Cactus for 
coupled spacetime and relativistic hydrodynamics evolution based on Riemann 
solvers. M. Miller implemented the key science module (MAHC) for GR3D, 
with further contributions from other members of the Grand Challenge H, |82| , 
see Sec. fj] for an application. Finally, Cactus is of course our platform for the 
binary black hole collisions reported on in Sec. ^| and the strong gravitational 
wave evolutions of Sec. ||. In this case, the code "BAM" originally developed in 
|33[ contributed the multigrid elliptic solver for the initial data and for 
maximal slicing, and the Mathematica scripts of BAM were used to generate C 
code for the BSSN evolution. For analysis, an apparent horizon finder imple- 
mented by M. Alcubierre was used (there is also one available by C. Gundlach, 
[|73"|), and the wave extraction routines by G. Allen A much larger number 
of individuals than is apparent from the above citations has contributed to "the" 
Cactus code, see J75|. At this moment the Cactus 3.2 CVS repository lists 88 
thorns, which range from private and under development to stable and public. 



5 Grazing collision of black holes 

The first crude but truly (3+l)-dimcnsional binary black hole simulation can 
be summarized as follows Q. The approach taken was to address each of 
the items listed in the skeleton for evolution problems of black holes of Sec. 



3.3| in the simplest possible manner that still allowed us to combine all the 
ingredients to a complete implementation. Initial data for two black holes, each 
with linear momentum and spin, is constructed using the puncture method ||66f 
(see also [84| ) , in which the internal asymptotically flat regions of the holes are 
compactified so that the numerical domain becomes R 3 . By construction the 
initial data is conformally flat. The evolution is performed with the original 
ADM equations and a leapfrog finite difference scheme. Maximal slicing and 
vanishing shift is chosen, i.e. physical singularities are avoided, while coordinate 
singularities typically do not occur for this elliptic slicing condition. At the 
outer boundary, the ADM variables are held constant, which works well for the 
achievable run times because a fixed mesh refinement of nested boxes (with finer 
resolution at the center) is used, and for large radii the lapse can approximate 
the Schwarzschild lapse for which Schwarzschild data would remain static. An 
important insight is that the puncture method, which can be made rigorous for 
initial data, can be numerically extended to the evolution equations so that no 
special inner boundary is present [ [4l| . Analysis is restricted to apparent horizon 
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finding with a curvature flow method. These methods allow one to evolve for 
about 7M, which is sufficient to observe the merger of the apparent horizons, 
but too short for wave extraction. 

Currently, binary black hole mergers are simulated by our AEI/NCSA/- 
WashU/Palma collaboration, and in this section I summarize some of the pre- 
liminary results. These simulations build on [^6|, Q and introduce various 
improvements. On the technical side, for high-performance collaborative com- 
puting the code is implemented with Cactus 3.2. An improved apparent horizon 
finder is now available Jt3[ . The comparatively slow maximal slicing can in many 
situations be replaced by "1+log" slicing (cmp. Sec. 3J2). No mesh refinement is 
used, but the outer boundary is treated with a radiative (Sommerfeld) boundary 
condition. The above plus the BSSN evolution system as given in Q with a 
3-step iterative Crank-Nicholson (ICN) scheme, allow run times of up to 30M 
for grazing collisions, compared to 7M for previous runs, and up to 50M for 
simpler data sets. 

The important new result is that now for the first time the extraction of wave 
forms becomes possible with the methods tested in Q . Let us discuss a concrete 
example. For initial data we choose the punctures of each hole on the y-axis 
at ±1.5, masses mi — 1.5 and m 2 — 1, linear momenta Pi 2 = (±2,0,0), and 
spins 5*1 = (—1/2, 1/2,0) and S 2 = (0, 1, 1) (all units normalized by m 2 = 1). 
The numerical grid has 385 3 points with grid spacing 0.2, which puts the outer 
boundary for a centered cube at a coordinate value of about 38. The initial 
ADM mass is M — 3.11, so the outer boundary is at about 12 M (solving 
the constraints for the "bare" parameters increases the mass over the Brill- 
Lindquist vanishing spins and momenta value of mi + m^). The total angular 
momentum is J — 6.7, which corresponds to an angular momentum parameter 
of a/M = J/M 2 = 0.70. 

The black holes start out with separate marginally trapped surfaces forming 
the apparent horizon (although it may well be that they have a common event 
horizon). Fig. ^ shows the formation of a single marginally trapped surface 
surrounding the initial inner marginally trapped surfaces. The apparent horizon 
is defined by a type of minimal surface equation (e.g. [[73)), and does not evolve 
continuously, rather a new "minimal" surface appears in a new location. The 
shading of the surfaces indicates the Gauss curvature on the surfaces. The area 
of the apparent horizon increases in these coordinates because grid points are 
falling into the black hole. In Fig. ||, two frames near the merger are shown 
together with isosurfaces of Re"04 as a wave indicator. 

From the evolution, we can obtain an energy balance by comparing the 
energy carried by the various modes of the gravitational waves |}5| with the 
difference in mass between the initial slice and the final black hole. For the 
latter, the apparent horizon mass is Mah = [(Mir) 2 + >/ 2 /(2M ir ) 2 ] 1 / 2 with 
(M ir ) 2 = A ah I (167r) and Aah the numerically determined area of the apparent 
horizon of the final black hole. During the evolution, Aah reaches a plateau, 
but then starts drifting upwards as the grid stretching becomes more severe. 
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Figure 2: The evolution of the apparent horizon during a grazing black hole 
collision (t = 9, 11, 16, 21). 




With the plateau value for A ah and Madm = M = 3.11 we obtain 



M ADM - M AH = 0.03 w O.OIMadm- (10) 

A rough estimate for the radiated energy in all modes until t = 30M, with the 
extraction radius rather close to the system at 8M, is 

M R ad = 0.007 - O.OOSMadm- (H) 

Even with all the current restrictions on accuracy coming from resolution, grid 
size, boundary treatment, and grid stretching, this energy balance can be con- 
sidered to be a first physics result for such a grazing collision. We learn that for 
this data set roughly 1% of the total energy is emitted in gravitational waves. 
Clearly, a thorough parameter space study of such configurations is of inter- 
est. To make contact with astrophysical situations, more realistic initial data is 
probably needed (which ideally would be derived from the slow inspiral of the 
two black holes). A detailed report is in preparation. 



6 Gravitational collapse of gravitational waves 



One way to probe general relativity in the highly non- linear regime, which should 
also share some of the strong wave features of the grazing collision, is certainly 
through the gravitational collapse of gravitational waves to a black hole. As 
briefly mentioned in Sec. 2.1, one scenario is that of critical collapse p3|. One 



can construct a one-parameter family of initial data, and examine the region 
near the "critical" value for that parameter at which a black hole does or does 
not form. Not much is known in 3+1 dimension s p4| , and the only study in 
axisymmetry is that of Abrahams and Evans p2| |25|] for gravitational waves. 

In this section I want to briefly discuss first results for non-axisymmetric 
collapse, cmp. [^9|. We take as initial data a pure Brill type gravitational wave 
[p5| , later studied by Eppley pB, pTj an d others [p8| . The metric takes the form 



ds 2 



tf 4 



■P-'i 



{dp 2 + dz 2 ) + p 2 



■>& 4 ds 



where q is a free function subject to certain boundary conditions. Following 
[33], jsijf , we choose q of the form 



(12) 
1 



a p 2 e 



1 



P 



(1 + P 2 ) 



cos 



(n</>) 



(13) 



= p 2 + z 2 and n is an 



where a,c are constants (a different from Sec. |^), r 
integer. For c = 0, these data sets reduce to the Holz |8£j] axisymmetric form, 
recently studied in three-dimensional Cartesian coordinates in preparation for 
the present work |fT3f . Taking this form for q, we impose the condition of 
time-symmetry, and solve the Hamiltonian constraint numerically in Cartesian 
coordinates. An initial data set is thus characterized only by the parameters 
(a, c, n). For the case (a, 0,0), we found in |73 that no apparent horizon exists 
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in initial data for a < 11.8, and we also studied the appearance of an apparent 
horizon for other values of c and n. 

For evolutions, we found that the BSSN system as given in EBj with maximal 
slicing, a 3-step ICN scheme, and a radiative boundary condition is sufficiently 
reliable even for the strong waves considered here. The key new extensions 
to previous BSSN results are that the stability can be extended to (i) strong, 
dynamical fields and (ii) maximal slicing, where the latter requires some care. 
Maximal slicing is defined by vanishing of the mean extrinsic curvature, K=0, 
and the BSSN formulation allowed us to cleanly implement this feature numer- 
ically in contrast with the standard ADM equations. 

As discussed in axisymmetric data with a — 4 is subcritical, that is 
the imploding part of the wave disperses again, leaving flat space in a non- 
trivially distorted coordinate system. An amplitude of a = 6 gives a super- 
critical evolution as indicated by the formation of an apparent horizon. The 
"cartoon" method |9(| to perform axisymmetric calculations in Cactus using 
three-dimensional Cartesian stencils on a two-dimensional slab allowed us to 
close in on the critical region near a = 4.6, but work on detection of critical 
phenomena is still in progress. 

Fig. U shows the development of the data set (a=6, c=0.2, n=l), which has 
reflection symmetry across coordinate planes. The initial ADM mass of this data 
set turns out to be Madm = 1-12. Fig. ^a shows a comparison of the apparent 
horizons of this three-dimensional and the previous axisymmetric cases at £=10 
on the x-z plane. The mass of the three-dimensional apparent horizon case is 
larger, weighing in at Mah=0-99 (compared to Mah{2D) = 0.87). 

In Fig. ^> we show the {Z=2,m=0} wave form of this three-dimensional 
case, compared to the previous axisymmetric case. The c = 0.2 wave form has 
a longer wave length at late times, consistent with the fact that a larger mass 
black hole is formed in the three-dimensional case. Figs, ^c and |]d show the 
same comparison for the {/=4,m=0} and {l=2,m—2} modes respectively. Notice 
that while the first two modes are of similar amplitude for both runs, the three- 
dimensional -{7=2,771=2} mode is completely different; as a non-axisymmetric 
contribution, it is absent in the axisymmetric run (in fact, it does not quite 
vanish due to numerical error, but it remains of order 10~ 6 ). We also show a 
fit to the corresponding quasi normal modes of a black hole of mass 1.0. The 
fit was performed in the time interval (10,36), and is noticeably worse if the 
fit is attempted to earlier times, showing that the lowest quasi normal modes 
dominate at around 10. The early parts of the wave forms t < 10 reflect the 
details of the initial data and BH formation process. This is especially clear in 
the {7=2, m=2} mode, which seems to provide the most information about the 
initial data and the three-dimensional black hole formation process. 

7 Minimal distortion shift 

As a final example for recent advances in numerical relativity simulations, let me 
mention shift conditions in (3+l)-dimensional relativity. The first preliminary 
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Horizon position at t = 1 Even mode, 1=2 m=0 
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Figure 4: a) The solid (dotted) line is the apparent horizon for the 3D data 
set (6, 0.2, 1) ((6, 0, 0)) at t=W on the x-z plane, b) The {Z=2,m=0} wave form 
for the 3D (6, 0.2, 1) case at r — 4 (solid line) is compared to axisymmetric 
(6,0,0) case (dotted line). The dashed line shows the fit of the 3D case to 
the corresponding mode for a black hole of mass 1.0. c) Same comparison 
for the {Z=4,777=0} wave form, d) Same comparison for the non-axisymmetric 
•{7=2,771=2} wave form. 
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test of a dynamically computed minimal distortion shift can be found in p0| 
for a Schwarzschild black hole on a 3+1 Cartesian grid, which is still the only 
example with black hole excision. Computational domains with holes pose a 
technical problem for the elliptic solver, which certainly will be solved (see for 
example once excision runs demand dynamic shifts. 

A non-vanishing shift plays an important role in calculations that involve 
orbiting black holes or neutron stars, e.g. in post-newtonian calculations or 
Newtonian hydrodynamics for neutron stars. The freedom in the shift vector 
can in principle be used to obtain corotating coordinates or partially corotating 
coordinates (to counter frame dragging). A variational principle to minimize 
coordinate shear leads to the minimal distortion family of shift conditions, see 
[ p2| . Introducing again a conformal factor such that the conformal metric g ab 
has unit determinant, one can minimize 

S[(3]= J\d t g\ 2 dV = J g ac g bd d t g ab d t9cd y^d 3 x, (14) 

which gives a vector elliptic equation for (3 a , 

(A^r = 2D b (a(K ab - g ab K/3)), (15) 
(A,/3) Q = D b D b (3 a + D b D a /3 b -~D a D b p b . (16) 

Note that if there exists a rotational Killing vector, minimal distortion can be 
trivially obtained 62J, hence such shift conditions begin playing a non-trivial 



role only when one moves beyond axisymmetric simulations (see also 1 93 , M , 
and in particular |95| for spacetimes with approximate Killing vector fields). 

There are now three examples for the application of dynamical shift con- 
ditions to binary neutron star simulations, which share the feature that with 
vanishing shift the code fails after far less than an orbit, while with minimal 
distortion shift for the first time fully relativistic simulations of one or more 
orbits become possible. In ]96"| , |97j |, minimal distortion is approximated in a 
way that decouples the three equations but maintains key features. Preliminary 
experiments have also been performed for the NASA Neutron Star Grand Chal- 
lenge using Cactus, the hydrodynamics module or "thorn" MAHC |8l], ^2|, and 
the author's implementation of the vector Laplace operator (|l^) in BAM. The 
full minimal distortion equations are solved. One choice of initial data is that 
of irrotational neutron star binaries provided by the Meudon group (|98|, |99|, 
polytropic equation of state with 7 = 2, K = 0.03c 2 /p nuc , M\ — M2 — 1.6M S0 ;, 
M/R = 0.14, d = 41fcm). Fig. || shows four frames of an evolution project that 
was implemented and carried out this summer by M. Miller, N. Stergioulas, and 
M. Tobias. Without shift, the simulation crashes before less than l/10th of 
an orbit is completed, with shift one observes about 3/4 of an orbit before the 
code fails when the two neutron stars merge. These first results can probably 
be improved significantly, but they already serve as a proof that non- vanishing 
shift is beneficial. 
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Figure 5: Binary neutron star orbit with minimal distortion shift. The neutron 
stars are represented by an isosurface at 1/10 of the central density, the shift 
vector by the arrows. 
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8 Conclusion 



It is perhaps surprising how little has been achieved to date by numerical sim- 
ulations of the Einstein equations for the two body problem. After all, the 
Einstein equations have been extensively studied for more than 80 years, and 
nowadays modern computational physics has successfully treated the partial 
differential equations of a large number of evolution problems. Why is it not 
possible to "simply solve" the problem with standard numerical methods on a 
big computer? To recall some of the issues raised above, (i) the Einstein equa- 
tions do not lead to a unique or preferred set of 3+1 evolution equations, with 
an automatically stable numerical implementation, (ii) choosing lapse and shift 
is intricately coupled to the evolution, (iii) black holes pose a special challenge 
due to their singularities. 

As a result, black hole simulations in numerical relativity still have to be 
called rather limited. Either special simplifications are introduced (axisymmc- 
try, null coordinates adapted to single black holes), or the achieved numeri- 
cal runtime is a limiting factor compared to the lowest quasi-normal ringing 
period of about 17M. (3+l)-dimensional black hole evolutions with singular- 
ity avoiding slicing last to about 30A/ for simple data sets starting from time 
symmetry (vanishing extrinsic curvature) |2^, ^0|. The first evolution of truly 
three-dimensional binary black hole data (two black holes with spin and linear 
momentum) was performed in 1997 [4l]| , crashing at 7M, which allowed tracking 
the merging of apparent horizons but not wave extraction. Considering that the 
first 3+1 simulations of Schwarzschild were reported in 1995, one can certainly 
call the recent simulations of Sec. ^| with wave extraction and a run time of 
about 30Af a significant step forward. 

Several methods are under intense investigation that should allow us to 
evolve for hundreds of M or even longer. Here we mentioned black hole ex- 
cision, improved evolution schemes, and shift conditions. Especially excision 
is expected to be essential. For the purpose of wave extraction, the schemes 
involving future null infinity are of particular interest. Furthermore, astrophys- 
ically more realistic initial data is needed as input for the above methods before 
we can make contact with gravitational wave astronomy. 

How close is numerical relativity to the accurate prediction of gravitational 
wave forms for binary events fllPPf ? The post-newtonian and the close- limit 
approximations are probably in good shape, but full numerical relativity will 
require two or more years to get ready. An introductory statement often heard 
during the last two decades is that one essential task of numerical relativity 
is to provide a catalog of wave forms which is essential for gravitational wave 
detection. This has changed. Numerical relativity will be essential in wave 
analysis, producing models for astrophysical scenarios that relate the wave forms 
to configuration parameters. For the detection as such, however, the task of 
producing a complete catalog appears to be too hard, and in particular, not a 
very sensible one. Note that matched filtering gives roughly a factor 5 in signal 



to noise for wave detection [101, [1021 . Recently, the advantage of the perfect 



catalog over the best "blind" numerical methods has been reduced to a factor 
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of 2, e.g. [ 103 1 . This still corresponds to a factor of 100 in observable event rate, 
but on the other hand optimal matched filtering is assumed in this estimate. 
The emphasis in numerical relativity should therefore be shifted more towards 
producing reliable statements about global features of mergers as opposed to 
detailed wave forms. Predicting the duration of mergers, total energy emission, 
frequency range and frequency distribution of the signal will be more useful to 
methods as described in [101, |10^_ 103 1 and also more attainable in the near 
future. The black hole runs of Sec. |5| are being performed with this goal in mind. 

It is a pleasure to thank E. Seidel and all the members of the numerical relativity 
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